%% Fig12a_Replication_Improved_Debug.m
% 复刻论文 Fig.12(a)：RM-DCSK 各混沌映射 BER vs. SNR 曲线
% 带有 Debug 步骤，便于检查中间环节是否正确。
% 注意：参数、帧结构、归一化处理均参考论文描述，若结果不理想，请检查每一步

clear; clc; close all;
rng(0);  % 固定随机种子，确保结果可重复

%% 调试开关：设置为 true 时输出调试信息和图形
debug = true;

%% 1. 参数设置
SNR_dB = 0:2:30;         % SNR 扫描范围（dB），可根据需要调整步长
N = 50;                  % Spread factor (每帧参考段或数据段长度)
numBits = 1e4;           % 发送比特数；调试时可先使用 1e4~1e5 测试

% 定义混沌映射名称（论文中 2D 和 1D 对比）
chaosMapNames = { '2D-MLM', '2D-MSM', '2D-MTM', '2D-SMM', ...
                  'Logistic', 'Sine', 'Tent' };

% 对应参数结构，注意参数与论文一致
paramsMap = { struct('mu',0.2, 'k',2.08, 'x0',0.2, 'q0',0.1), ... % 2D-MLM
              struct('mu',-0.1, 'k',1.88, 'x0',0.8, 'q0',0.6), ...  % 2D-MSM
              struct('mu',-0.6, 'k',1.78, 'x0',0.8, 'q0',-0.8), ...  % 2D-MTM
              struct('mu',0.5, 'k',2.33, 'x0',1, 'q0',1), ...        % 2D-SMM
              struct('mu',3.68, 'x0',0.2), ...                       % Logistic (1D)
              struct('mu',0.93, 'x0',0.2), ...                       % Sine (1D)
              struct('mu',1.49, 'x0',0.2) };                         % Tent (1D)

% RM-DCSK 帧结构：每帧长度 = 2*N（前 N 参考段，后 N 数据段）
% 为保证各映射在同一 SNR 下起始点一致，对混沌序列做统一归一化

%% 2. 主仿真循环
numMaps = length(chaosMapNames);
BER_matrix = zeros(numMaps, length(SNR_dB));

for mapIdx = 1:numMaps
    mapName = chaosMapNames{mapIdx};
    fprintf('Simulating map: %s\n', mapName);
    % 获取对应参数结构
    params = paramsMap{mapIdx};
    
    % 产生混沌序列：每个比特需 2*N 个采样，总长度 = 2*N*numBits
    chaosSeq = generateChaosSequence(mapName, numBits, N, params);
    
    % Debug: 绘制混沌序列前 1000 点并输出均值、标准差
    if debug
        figure('Name', ['Chaos Sequence: ', mapName]);
        plot(chaosSeq(1:min(1000, end)));
        title(['前 1000 个采样: ', mapName]);
        xlabel('样本点'); ylabel('幅值');
        fprintf('Map: %s, chaosSeq 均值=%.4f, 标准差=%.4f\n', ...
                mapName, mean(chaosSeq), std(chaosSeq));
    end
    
    % 对混沌序列进行归一化（统一平均功率），确保均值为0，标准差为1
    chaosSeq = (chaosSeq - mean(chaosSeq)) / std(chaosSeq);
    
    % Debug: 输出归一化后信号的均值和标准差
    if debug
        fprintf('After normalization: mean=%.4f, std=%.4f\n', mean(chaosSeq), std(chaosSeq));
    end
    
    % 生成随机比特，取值 {-1, +1}
    dataBits = randi([0,1], 1, numBits)*2 - 1;
    
    % RM-DCSK 编码：构建 txSignal
    txSignal = RM_DCSK_Encode(dataBits, chaosSeq, N);
    
    % Debug: 检查 txSignal 的均方功率
    txPower = mean(txSignal.^2);
    fprintf('Map: %s, txSignal 均方功率=%.4e\n', mapName, txPower);
    
    % 对每个 SNR 进行仿真
    for iSNR = 1:length(SNR_dB)
        snrVal = SNR_dB(iSNR);
        
        % 通过 AWGN 信道 (使用 'measured' 模式)
        rxSignal = awgn(txSignal, snrVal, 'measured');
        
        % RM-DCSK 解码：得到接收比特 rxBits
        rxBits = RM_DCSK_Decode(rxSignal, N);
        
        % 计算 BER
        BER = sum(rxBits ~= dataBits) / numBits;
        BER_matrix(mapIdx, iSNR) = BER;
        
        % Debug: 可打印前几个 SNR 下的 BER 观察变化
        if debug && iSNR <= 2
            fprintf('Map: %s, SNR=%.1f dB, BER=%.4e\n', mapName, snrVal, BER);
        end
    end
end

%% 3. 绘图：BER vs. SNR (半对数坐标)
figure('Name','Fig.12(a) Replication - Improved Debug');
hold on; grid on;
markers = {'o','s','^','d','+','x','*'};  % 定义绘图标记

for mapIdx = 1:numMaps
    semilogy(SNR_dB, BER_matrix(mapIdx,:), 'LineWidth',1.5, ...
        'Marker', markers{mapIdx}, 'MarkerSize',6, ...
        'DisplayName', chaosMapNames{mapIdx});
end
xlabel('SNR in dB');
ylabel('BER');
title('信噪比-误码率曲线');
legend('Location','Best');

%% ================== 子函数 ==================

%% 产生混沌序列
function chaosSeq = generateChaosSequence(mapName, numBits, N, params)
    % 每个比特帧长度 = 2*N, 总点数 = 2*N*numBits
    totalLen = 2 * N * numBits;
    chaosSeq = zeros(1, totalLen);
    
    switch mapName
        case '2D-MLM'
            mu = params.mu; k = params.k; x = params.x0; q = params.q0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                [x, q] = iteration_2DMLM(x, q, mu, k);
            end
        case '2D-MSM'
            mu = params.mu; k = params.k; x = params.x0; q = params.q0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                [x, q] = iteration_2DMSM(x, q, mu, k);
            end
        case '2D-MTM'
            mu = params.mu; k = params.k; x = params.x0; q = params.q0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                [x, q] = iteration_2DMTM(x, q, mu, k);
            end
        case '2D-SMM'
            mu = params.mu; k = params.k; x = params.x0; q = params.q0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                [x, q] = iteration_2DSMM(x, q, mu, k);
            end
        case 'Logistic'
            mu = params.mu; x = params.x0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                x = mu * x * (1 - x);
            end
        case 'Sine'
            mu = params.mu; x = params.x0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                x = mu * sin(2*pi*x);
            end
        case 'Tent'
            mu = params.mu; x = params.x0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                if x < 0.5
                    x = mu * x;
                else
                    x = mu * (1 - x);
                end
            end
        otherwise
            error('Unknown chaos map name!');
    end
end

%% RM-DCSK 编码
function txSignal = RM_DCSK_Encode(dataBits, chaosSeq, N)
    % 每比特对应一帧，帧长 = 2*N
    L = length(dataBits);
    txSignal = zeros(1, 2*N*L);
    idx = 1;
    for i = 1:L
        % 前 N 点为参考段
        ref = chaosSeq(idx : idx+N-1);
        % 后 N 点为数据段，乘以比特值
        dat = dataBits(i) * chaosSeq(idx+N : idx+2*N-1);
        txSignal((2*N*(i-1)+1):(2*N*i)) = [ref, dat];
        idx = idx + 2*N;
    end
end

%% RM-DCSK 解码
function rxBits = RM_DCSK_Decode(rxSignal, N)
    % 每帧 2*N 点，采用相关器判决
    L = length(rxSignal) / (2*N);
    rxBits = zeros(1, L);
    idx = 1;
    for i = 1:L
        ref = rxSignal(idx : idx+N-1);
        dat = rxSignal(idx+N : idx+2*N-1);
        correl = sum(ref .* dat);
        if correl > 0
            rxBits(i) = +1;
        else
            rxBits(i) = -1;
        end
        idx = idx + 2*N;
    end
end

%% 迭代方程：2D-MLM
function [x_next, q_next] = iteration_2DMLM(x, q, mu, k)
    x_next = mu * x * (1 - x) + k * cos(q) * x;
    q_next = q + x;
end

%% 迭代方程：2D-MSM
function [x_next, q_next] = iteration_2DMSM(x, q, mu, k)
    x_next = mu * sin(2*pi*x) + k * cos(q) * x;
    q_next = q + x;
end

%% 迭代方程：2D-MTM
function [x_next, q_next] = iteration_2DMTM(x, q, mu, k)
    if x < 0.5
        x_next = mu * x + k * cos(q) * x;
    else
        x_next = mu * (1 - x) + k * cos(q) * x;
    end
    q_next = q + x;
end

%% 迭代方程：2D-SMM
function [x_next, q_next] = iteration_2DSMM(x, q, mu, k)
    x_next = mu * x + k * cos(q) * x;
    q_next = q + x;
end
